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Abstract 



Aims. In this paper we study density cusps that may contain central black holes. The actual co-eval self-similar growth would not 
distinguish between the central object and the surroundings. 

Methods. To study the environment of a growing black hole we seek descriptions of steady 'cusps' that may contain a black hole and 
that retain at least a memory of self-similarity. We refer to the environment in brief as the 'bulge' and on smaller scales, the 'halo'. 
Results. We find simple descriptions of the simulations of collisionless matter by comparing predicted densities, velocity dispersions 
and distribution functions with the simulations. In some cases central point masses may be included by iteration. We emphasize 
that the co-eval self-similar growth allows an explanation of the black hole bulge mass correlation between approximately similar 
collisionless systems. 

Conclusions. We have derived our results from first principles assuming adiabatic self-similarity and either self-similar virialisation 
or normal steady virialisation. We conclude that distribution functions that retain a memory of self-similar evolution provide an 
understanding of collisionless systems. The implied energy relaxation of the collisionless matter is due to the time dependence. Phase 
mixing relaxation may be enhanced by clump-clump interactions. 

Key words, theory-dark matter-galaxies:haloes-galaxies:nuclei-black hole physics-gravitation. 



1. Introduction 

In a previous work (Henriksen et al., referred to here as paper 
|I] and references hereafter) we discussed the relation between 
the formation of black holes (hereafter BH) and of galaxies. We 
presented distribution functions (DF) for the co-eval formation 
of spherically symmetric cusps and bulges, as formed by radially 
infalling, collisionless matter (e.g. dark matter or stars). We also 
discussed the influence of a centrally dominant mass, including 
a point mass BH. 

Observations ( |Kormendy & Richstone 1995 



Mag orrian et al. 1998| IFerrase & Merritt 2000 

Gebhardt et al. 2000) have established a strong correlation 
between the BH mass and the surrounding stellar bulge mass (or 
velocity dispersion), which we take to be an indication of co- 
eval growth. Such growth occurs in part during the dissipative 
baryon accretion by BH 'seeds' in the AGN (Active Galactic 
Nuclei) phase, but there is as yet no generally accepted scenario 
for the origin of the seeds. Moreover, recent suggestions of very 
early very supermassive BHs (e.g. IKurk et al. 20071 1. together 
with changes in the normalization of the BH mass-bulge mass 
proportionality (relatively larger black holes at high red shift) 
(e.g. IMaiolino et al. 20071 ). suggest an alternate early growth 
mechanism. Recently |Peirani & de Freitas Pacheo (2 008) have 
studied the possible size of the dark matter component in BH 
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masses. They deduce that between 1% and 10% of the black 
hole mass could be due to dark matter. 

We explored this possibility in the case of spherically sym- 
metric radial infall (paper I) and in this paper we will extend it 
to anisotropic infall. We use the same technique of inferring rea- 
sonable distribution functions for collisionless matter from lim- 
its of the time dependent Collisionless Boltzmann (CBE) and 
Poisson set, and are aided by some high resolution simulations 
of the evolution without a dominant central mass. We add a 
dominant central mass either analytically or by iteration. Just 
as in the radial case, the loss cones are not empty for substan- 
tial growth. The relaxation of collisionless matter is due to the 
temporal evolution (including the radial orbit instability) and, 
in addition, to possible 'clump-clump' (two clump) interactions 
dHenriksen 20091 IMacMillan & Henriksen 20031) . 

We use, in this paper as in the previous paper ©, the 
Carter-Henriksen (Carter & Henriksen 19911) procedure to ob- 
tain a quasi-self-similar system of coordinates (Henriksen 
2006a,2006b, hereafter lH20Tj6a1 |H2006bl This allows writing 
the CBE-Poisson set with explicit reference to possible transient 
self-similar dynamical relaxation. In this way we can remain 
'close' to self-similarity just as the numerical simulations appear 
to do. 

Studies of BH-density cusps originated with the problem of 
BH feeding dPeebles 1972IIBahcall & Wolf 19761), and with the 
notion of adiabatic growth (|Young 1980| IQuinlan et al„ 1995] 
MacMillan & Henriksen 2002][ Observations of the cen- 
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tral Milky Way have detected a mainly isotropic den- 
sity cusp with logarithmic power in the range —1.1 ±0.3 
dGillessen et al. 2009j l. This is flatter than the adiabatic limit 
(MacMillan & Henriksen 2002) and in any case the adiabatic 
growth scenario does not produce the BH bulge correlations 
(ibid). All this has spurred the investigation of co-eval dy- 
namical growth instead of adiabatic growth. Central cusps 
flatter than —1.5 can be created by tight binary BH systems 
formed in mergers that 'scour' the stellar environment (e.g. 
IMerritt & Szell 20061 INakano & Makino 1999t . This process 
can produce log slopes as flat as —1 or even —0.5, and it is 
supported by a strong correlation between nuclear BH mass and 



central luminosity deficit dKormendy & Bender 20091. 

However the correlation in itself only implicates the influ- 
ence of the BH. It does not necessarily require the merger his- 
tory, which in any case is unlikely to be the same for different 
galaxies. Consequently we explore in this paper whether cusps 
as flat as those resulting from scouring might also be produced 
during the dynamical formation of the black hole. 

We begin the next section with a summary of the general for- 
mulation in spherical symmetry. Subsequently we study a sys- 
tem comprised of anisotropic non-radial orbits in spherical sym- 
metry. Finally we give our conclusions. 

2. Dynamical equations and previous results 

We will use the formulation of H2006a, wherein we transform 
to infall variables the collisionless Boltzmann and Poisson equa- 
tions for a spherically symmetric anisotropic system. We begin 
with the 'Fujiwara' form (e.g. Fujiwara 1983 ), namely 



^/^ df 
dt +Vr dr 



j 2 d<P\ df 



= 0, 



d_ 



,<94> 



r ' dr ) dv r 
= An 2 G I f{r,v r J 2 )dv r dj 2 . 



(1) 

(2) 



Here / is the phase-space mass density, <t> is the 'mean' field 
gravitational potential, j 2 is the square of the specific angular 
momentum and other notation is more or less standard. 

The transformation to infall variables takes the form (e.g. 
H2006i 



R = re- aT '\ 
;2„-(4/a-2)aT 



Y=v r e-^ a -^ aT , 
e aT = at, 



P(R,Y,Z;T) = e^ a -^ aT nf(r,v r ,j 2 ;t). 
y{R-T)=e- 2{l l a -^ aT <S>{r), &{R;T) = p(r,t) e - 2aT . 



(3) 



The passage to the self-similar limit requires taking dj = 
when acting on the transformed variables. Thus the self-similar 
limit is a stationary system in these variables, which is a state that 
we refer to as 'self-similar virialisation' (Henriksen & Widrow 
1999, hereafter |HW 19991 ILe Delliou 200TT l. The virial ratio 
2K/\W\ is a constant in this state (although greater than one; 
K is kinetic energy and W is potential), but the system is not 
steady in physical variables as infall continues. 

The single quantity a is the constant that determines the dy- 
namical similarity, called the self-similar index. It is composed 
of two separate reciprocal scalings, a in time and 8 in space, in 
the form a = a/ 8. As it varies it reflects all dominant physical 
constants with dimensions of mass as well as of length and time. 



This is because the mass reciprocal scaling ji has been reduced 
to 35 — 2a in order to maintain Newton's constant G scale in- 
variant (e.g. lH2006al ). 

We assume that time, radius, velocity and density are mea- 
sured in fiducial units r g /v ,r , v and p r) respectively. The unit 
of the DF is f and that of the potential is v 2 . We remove con- 
stants from the transformed equations by taking 



fo=Po/v 3 , v 2 =\%Gp r 2 . 



(4) 



These transformations convert equations (H},© to the re- 
spective forms 



-d T P-(3/a-l)P+(---)d R P 
a a a 

-f(l/a-l)Y + U^-^))dyP-(4/a-2)Zd z P = 



and 



1 d (r 2 ^ 1 O 



This integro-differential system is closed by 



1 



PdY dZ. 



(5) 



(6) 



(7) 



This completes the formalism that we will use to obtain the 
results below. The 'cusps' we describe there will generally end 
in what is the central 'bulge' surrounding the black hole, rather 
than in the black hole itself. 



3. Spherically symmetric steady anisotropic bulges 
and cusps 

We follow the same pattern of discussion in this section that we 
used previously (paper 1]) for radial orbits. We begin in this sec- 
tion with steady bulge-black hole systems that retain a memory 
of the prior nearly self-similar relaxation. The relaxation pro- 
ceeds by way of the relation 



dE 
~dt 



d<P 



(8) 



(with the appropriate total energy and potential energy). This in- 
cludes clump-clump interactions and the radial orbit instability 
and produces finally the coarse grained (i.e. steady H2006aJ) sys- 
tem. 

The steady self-similar cusps were found in (Henriksen & 
Widrow 1995, hereafter HW95 ) except for the case where a = 1. 
However we can recover them here by using the same procedure 
that we used for radial orbits (see paper I ). We employ the char- 
acteristics of ©, together with the identity ^ = ^ + ^^^P. 
A combination of the characteristics leads to an expression for 
the scaled energy on a characteristic, namely 



dS 
ds 



— = -2(l/a-l)#+2(l/a- 



R dW dl 1 

a oR ds 



where 



Z 



(10) 
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We impose the steady state in equation ((9]) by setting 
d^/ds = and requiring as in paper I that *P oc RP with p = 
2(1 - a) so that 



Rd*¥ 
a~dk~ 



■2(l/a- 1)¥ = 0. 



(11) 



Consequently the characteristics of equation (|5j are seen to im- 
ply 



§ = (3/«-l)* 
£ = -2(l/a-l)*, 
^ = -(4/a-2)Z. 



(12) 



These combine to give the general self-similar steady DF in 
the form (a ^ 1 and a 7^ 3 /2) 



P = P(Z 0) k)^, 



where 



3 — a 
2(a-l)' 



(13) 



(14) 



The quantity Z c is the value of Z at s = aT = 0, that is j . 
However in the ab-initio steady self-similar analysis, this depen- 
dence does not appear (HW95 ). It is however present in general 
in this limit from time dependent phase, and we use it in the next 
paper in this series (paper III), on non-self-similar cusps. 

The physical form of the self-similar DF becomes, on using 
the scaling relations ©, 



nf = P{k)\E\ 



(15) 



where E = v 2 /2 + j 2 /(2r 2 ) 



<t>. This result has also been shown 
(e.g. IH2006al) to be the zeroth order DF in the coarse graining 
approach that becomes the exact steady state in the long time 
limit (a — > °°). It was also used by Kulessa and Lynden-Bell 
(Kulessa & Lynden-Bell 1992) in their study of the mass of the 
Galaxy. 

An earlier important paper by Stiavelli and Bertin 
dStiavelli & Bertin, 1985| l also proposed equilibrium distribution 
functions for elliptical galaxies. These were based in general on 
three integral models, but were reduced to our two familiar inte- 
grals in the case of spherical symmetry. In our case the motiva- 
tion for the equilibrium DF is different, because it derives from 
the self-similar dynamical growth. There are nevertheless some 
significant similarities, particularly as modified to an inverted 
energy distribution by Merritt, Tremaine and Johnstone (1989, 
hereafter lMTJI ). We discuss these points further below and in the 
conclusions to this paper. 

The density profile that accompanies this DF for a particular 
choice of P(k) can be shown by direct integration over phase 
space to be p °= r~ 2a , and subsequently the consistent potential 
is oc r 2 ^- a ) so long as a 7^ 1,3/2. 

The case a = 3/2 is excluded simply because it represents a 
point mass in a massless halo. However this is of interest pre- 
cisely when describing the environment of a dominant central 
mass, whether this be the halo of a central black hole or the halo 
of a coarse-grained collisionless bulge of stars. 



The general form of the DF becomes from equation (15[ 

nf = P(\E\j 2 )\E\ 3 ' 2 . (16) 

This gives a number density of massless particles « r~ 3 for a 
potential <t> = —M*/r. The rms radial velocity (the same as the 
radial velocity dispersion since v7 = 0) will generally be oc <J> for 
any power law choice of P(\E\ j 2 ). 

One can come close to imitating the Bahcall and Wolf 
(Bahcall & Wolf 19761 1 zero flux solution for a black hole cusp 
by choosing P(k) = F(\E\j 2 )/(\E\j 2 ) 5 / 4 . This yields the DF 



7Cf = F(\E\f 



\E\^_ 

'(,•2)5/4' 



(17) 



where F is an arbitrary function. In fact, they use for argument 
of their arbitrary function A = (1 — <? 2 )/2, where e is the orbital 
eccentricity, and for Keplerian orbits around a point mass M. 
this becomes \E\j 2 /(GM,) 2 . In our units GM. = 1. 

If this were exactly the Bahcall- Wolf solution it would be re- 
markable, since their solution is the result of collisional diffusion 
into the loss cone. Here however we are obliged to have the addi- 
tional dependence ( j 2 )~ 5 / 4 , which implies a divergent loss cone. 
Thus the necessary diffusion is simply assumed in this example 
in order to duplicate the energy and eccentricity dependence (ar- 
bitrary because of spherical symmetry). 

One can obtain the Bahcall- Wolf energy dependence from 
the general form (TT3T > only by setting a = 7/3, which is quite 
unrealistic and gives the wrong arbitrary dependence. We shall 
see in paper [III] that it is possible to choose a simple non-self- 
similar DF that does fall between these two cases and imitates 
the Bahcall- Wolf result more precisely. 

There are two simple limits of the form ([T81 for which there 
is some evidence in the numerical simulations of isolated dark 
matter bulges (MacMill an 20061 1. Just as in the radial case the 
following comparison with simulations suggest that this family 
of DF's may actually be realized. These limits are found by tak- 
ing first P to be a constant K, and then by taking it to be propor- 
tional to KT q . This yields respectively 



nf = K\E\", 
nf K 



(18) 
(19) 



where we have defined w = (3 — a) /{A — 2d). 

We expect (see discussion to follow) that the DF ( TT~8b will 
apply near the centre of a system where a < 1 . The self similar 
potential is then increasing with radius and so must be taken 
positive to obtain an attractive force. This allows us to calculate 
the corresponding density explicitly as 



p=4V2K(<t>)^B(\q\ -3/2,3/2). 



(20) 



Here B(x,y) is the standard Beta function, which may be ex- 
pressed in terms of the gamma function F(x) as T(x)T(y) /T(x + 

y). 

The Poisson equation shows that any power law solution for 
the potential with this density goes as <£> oc r 2 ( 1 -") and hence 



also p 



„-2a 



, both as they must for self-similarity. The only 



reason for giving this explicit form here is that it may be used in 
an iterative calculation of the transition from halo to black hole 
dominance. In such a calculation a point mass potential plus the 
halo potential is inserted to give a new density, which in turn 
may be used in the Poisson equation to give a new potential and 
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so on. The halo potential must dominate however so as to keep 
the net potential positive. 

As an example we assume that the potential has the form 
<t> = — C\/r + C2?" 2 ( 1 ~ a ) (Ci and C2 arbitrary positive constants) 
so that by d20j the density becomes p oc (C 2 - Ci/r 3 ~ 2a )r~ 2a . 
From the Poisson equation the new potential now follows as 



<p< 



2(l-a)(3-2fl) 



•(Cj-Ci In r)/r, 



(21) 



where C3 should be positive to maintain the point mass contribu- 
tion. The inner boundary of this expression is at the radius where 
it equals zero and the outer boundary would be inside the NFW 
(Navarro et al. 1996, hereafter INFW199 6) scale radius (which 
we may set equal to our scale radius r ). The radial mean square 
velocity v 2 (also the radial squared dispersion) is proportional 
to <t> and so we predict a minimum value just outside the radius 
r m where the central mass becomes dominant (i.e., <&(r m ) = 0), 
followed by a steady rise that tends to r 2 ( 1_a ). 

The density corresponding to this potential follows from (|20T > 



as 



(22) 



so that it has for a < 1 a rapid rise near r m followed by the self- 
similar regime wherein p °= r~ 2fl . The radius r m would be be- 
tween 1 and 10 pc in the halos of typical galaxies. 

Of course this iteration has not fully converged but subse- 
quent cycles will produce higher order terms. The result (f2Tb 
seems to indicate that the black hole mass is enhanced by an r~ 3 
cusp that peaks just outside r m . 

The radial velocity dispersion is generally equal to the root 
mean square radial velocity for the DF ( TT~8T > since the mean radial 
velocity is zero. This is proportional to 4> °= r 2 ' 1- "'. The squared 
tangential velocity dispersion however becomes 



erf =4> 



8B(| 9 |-5/2,5/2) n 2 B(\q\-2,2) 2 



3 B(\q\- 3/2.3/2) 8 -3/2,3/2)2 



(23) 



where q < —2 or 1 > a > 1 /3. The factor F(a) in this equation 
that multiplies \/<t> is shown in figure ([TJ as a function of a. 

Unlike the radial velocity dispersion, which should steadily 
increase with radius proportionally to <t> (which flattens as a — > 
1), the tangential velocity dispersion is expected to show a more 
precipitous drop as r tends towards the scale radius and a — >• 1 
as F(a) declines. The variation of a is only adiabatic with radius 
however, varying roughly as r 2 . We shall see evidence for this 
below. 

We turn first to an examination of the other limit ( fT9l . We 
expect this to apply in the outer region to which much angular 
momentum has been transferred by the bar due to the radial orbit 
instability (MacMillan et al. 2006, hereafter lMWH 20061 ). Hence 
we choose the regime where a > 1 and the energy/potential is 
negative. Iteration is not really relevant in this limit since we are 
far from the black hole, but because of its importance in calcu- 
lating mean quantities we give the density as 



?(3/2-w) 

p = — KI { W )r- w \<S>\^-»\ 

3/2 — w 



(24) 



Once again the Poisson equation has the consistent solution <t> < 
_ r 2 (i-«) an( j p ^ r -2« T ne integral Iq(w) is given by 



(25) 



0.75 — 




0.25 



Figure 1. The curve shows the dependence of the factor that 
multiplies <t> in the tangential velocity dispersion as a function of 
a. The adiabatic variation of a with r is approximately a = r 2 
(IH20071 >. where r is measured in units of the NFW scale radius 



where the lower limit must be greater than zero for convergence 
since w > 1 for a > 1 . 

The mean square radial velocity becomes in this limit 



2(3/2 -w) h 
5/2- w / 



(26) 



where, since a > 1, this decreases with radius. The tangential ve- 
locity dispersion has once again a more complicated expression 
namely 



< al >= |4>| 



2(3/2 -w)I (w- 1) 
5/2 — w h{w) 

/ 2 5 / 2 (3/2-w)/o(w-l/2) ^ 
^ (2-w)/ (w) 



(27) 



Although w(a) is slowly varying for 1 < a < 3/2, the factor mul- 
tiplying <t> declines to zero asiv-> 3/2. 

It is time however to compare these DF's chosen for their 
self-similar memory to the results of simulations, in order to jus- 
tify our attention. 

In references (112006a) and (Henriksen 2007, hereafter 
IH20071 ) it was concluded that a w 0.72 near the outer boundary 
of the relaxed region. Moreover we know from those papers, and 
more generally from extensive cosmological simulations, that 
a « 0.5 in the interior, well relaxed, region. This is referred to as 
adiabatic self-similarity (H2006a) since the index a dynamically 
evolves, but relatively slowly (a ^ r aF , a F = O(0.2), IH2007I) 
as already employed above. To match the simulation results we 
require a ss 3/2 in the vicinity of the NFW scale radius, after 
which there may be a tidal truncation to r~ 4 (Henriksen 2004, 
hereafter H2004ll 
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There is evidence in the simulations for the persistence of 
self-similarity, even in the case that most closely resembles an 
isolated cosmological halo (i.e. cosmological perturbations are 
included in the initial conditions). Thus in the cosmological-like 
halo simulation of (MacMilla n 20061 1. the virial mass and virial 
radius are found to grow in a power law manner according to 



t and f (we do not quote numerical errors, which are at 
the level of a few percent) respectively. The logarithmic density 
slope in the outer relaxed region is approximately — 1 .4 while 
the pseudo-density logarithmic slope is —2.16. 

The numerical predictions follow by using a = 0.72 from 
dH2007l l. The self-similar mass growth inside any growing ra- 
dius (fixed R) is (1HW 1999b °= f( 3 A<-2) oc f 216 , while r °= t l l a °c 
f 1 - 38 (see equation (0). The logarithmic density slope is pre- 
dicted to be —2a = 1.44, and the pseudo-density logarithmic 
power (e.g IH2007I Hansen 20041 is given by a - 3 = -2.28. 
These reasonable correspondences with the numerical simula- 
tion (only the predicted pseudo-density differs significantly from 
the simulated value) encourage us to adopt a DF as one of the 
above steady forms. The self-similar memory is incorporated ei- 
ther in a = 0.72 or in a = 0.5 in the relaxed region, which bound- 
ary is approximately (in fact the region where a < 1 extends 
beyond the scale radius by about a factor 2 in the simulation) 
coincident with the NFW scale radius. 

The velocity dispersion in the relaxed region is based on the 
DF ( TT8l and has been given in equation (123b for the tangential 
case, while the radial dispersion is proportional to \/<t> °< r l ~". 
Hence, away from the transition to the black hole region, we 
can test our predicted dispersions against the simulation . Using 
the value a = 0.72 this gives the radial (a>) velocity dispersion 
going as r° 28 which should be compared with figure (fJJ. The 
cosmological case is the right hand panel, which the phase space 
diagram (not shown) shows to be more relaxed than the initially 
unperturbed simulation on the left. The r 028 rise is a reasonable 
description of the numerical behaviour, especially as some adi- 
abatic evolution in a towards unity is to be expected. This will 
flatten the rate of rise Between 4 and 60 kpc, the radial disper- 
sion rises by about a factor 2 while the predicted value using 
a = 0.72 is 2.13. Allowing a to increase only to 0.75 on average, 
improves the fit. 

Eventually at large enough radius, according to the complete 
cosmological simulations, we require a 3/2 and so we predict 
a decline in the radial dispersion proportional to r~ x l 2 in this 
region. Figure (f2| indicates that this begins at about twice the 
indicated scale radius. Before this point but still outside the scale 
radius the decline is slower. In fact the simulation gives a « 5/4 
in this region so that the expected decline would be only as r~ '/ 4 . 

The tangential dispersion velocity has the same behaviour at 
small r but rolls over more rapidly. This is indicated very clearly 
in the graphs of the anisotropy parameter /3 = 1 — c 2 /a 2 on the 
panel below the velocity dispersion. This is consistent qualita- 
tively with the action of the factor F(a) defined in equation ( f23l . 
and displayed in figure ([TJ over the relevant range of a. The weak 
adiabatic dependence of a on r is not weak enough however to 
cover the range in r over which ax stays flat. There is some un- 
certainty in this factor and for a « r° 5 one arrives at nearly a 
factor two in radius as a varies from 0.72 to 0.8. However figure 
(Q]) shows that F has dropped by more than a factor two in this 
range, while V$ increases by only a factor of 1 .2 at best. The 
implied factor two decline is not seen. The only explanation is 
that a does not vary significantly in this region, so that there is 
little adiabatic evolution. 




r (kpc) 

Figure 2. The figure shows two distinct simulations from 
dMacMillan 20061 1. The right hand panels display the result (with 
continuing infall) of an isolated halo simulation starting from a 
set of particles perturbed by a cosmological spectrum of density 
and velocity fluctuations. The panel on the left indicates the same 
stage as developed from a non-perturbed halo. The upper line on 
the upper panel in each case is for a similar halo, but allowing 
only purely radial collapse. The middle line is the radial velocity 
dispersion while the lower line is the tangential dispersion. The 
lower panels display the anisotropy parameter for non-radial col- 
lapse. 



At the end of this plateau we have entered the region where 
a « 5/4 > 1 . The decline in v<& is however only as r~ 025 there, 
whereas the observed decline in figure (0 is more like r~ ! . Once 
again the difference between this behaviour and that of the radial 
dispersion must lie in the factor multiplying <t>. 

So we have found above that the implications of equations 
(fT8l and ( fT9] l are reasonable, but do we have any evidence con- 
cerning the DF itself? 

In the outer part of the relaxed region MacMillan (ibid, and 
figure O finds that the DF is mainly dependent on angular mo- 
mentum. This is evident from the figure by noting that numeri- 
cally g(EJ 2 )°< dM/dE in the range —100 to —200 energy units 
(MacMillan 2006). Hence from the definition 



f(E,f) 



1 



d l M 



g(E,f)dEdf 



(28) 



we infer that /°c d\ndM/dE/d j 2 . This is independent of energy 
in the outer region. Note that in terms of j 2 the figure is simply 
stretched by a factor of 2 in the j direction. 

We therefore fit the self-similar DF ( fT9l ) with a = 1.25 
to obtain / « (j 2 ) . This is in fact a reasonable fit to the 
DF in the less tightly bound region with angular momentum 
greater than about 500 units (figure [3}. Such a cut-off in an- 
gular momentum was contained in the Stiavelli and Bertin DF 
dStiavelli & Bertin, 1985| l where it was exponential. MTJ gave 
a heuristic justification for a cut-off in j and also wondered 
whether an exponential or a power-law cut-off was superior. 
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They also asked whether such behaviour arose naturally in the 
course of the dynamical formation of galaxies. Here our answer 
is in favour of a power-law cut-off following self-similar infall. 

An intermediate region where the logarithmic density slope 
varies from just above 2 to about 3 does occur in all simula- 
tions of collisionless matter. Moreover this region tends to co- 
incide with the passage from isotropy to radial orbits (e.g. |2}. 
According to the preceding remarks (i.e. the DF ( fT9T > with a > 1) 
this region can be characterized by roughly constant numbers of 
particles over a broad range of energy, but with numbers rapidly 
increasing towards a minimum in angular momentum. Some par- 
ticles would have too much angular momentum to enter the cen- 
tral region, while those below a critical angular momentum can 
populate the central regions. 

In the simulations this redistribution of angular momentum 
is attributed to the onset of the Radial Orbit Instability (ROI) 
in dMWH 20061) and the consequent formation of a bar. In the 
current spherically symmetric analysis, we are forced to patch 
these different regions together 'by hand', guided only by the 
different possibilities for a. 

In the inner relaxed region we expect the DF ( fTST l to apply 
with a ss 0.72. This gives an isotropic / E~ 4 . The simulations 
dNFW 19961 1 suggest that a — > 1/2 in the very centre of the sys- 
tem, which would give /pcE -V 2 . This is in accord with previ- 
ous discussions (e.g. H2006a). The fact that our energy is posi- 
tive, while the simulation energy is negative, reflects the differ- 
ent zero points for the potential. We use the centre of the system 
as the reference zero in our analysis, while the simulation uses 
infinity. 

Numerically, g(E) °= lEl -2,5 in this region so that 
d 2 M/dEdf should vary as E~ 6 - 5 (E~ 5 if a = 0.5). There is 
very little dependence on j 2 here so this prediction might be 
compared to the behaviour of dM/dE. This quantity is falling 
steeply beyond E = —300 (O, but the evidence is weak at the 
limit of resolution. 

At slightly higher angular momentum there is a rising de- 
pendence on j. This can be imitated with a self-similar solu- 
tion by taking P{k) oc \/k in equation (1151 1. which gives the DF 

/oc£-5(/)l/5. 

The above prolonged discussion suggests that a DF with a 
self-similar memory can be used to parametrise the phase space 
of these collisionless simulations. At least this is true if we ac- 
cept the adiabatic variability of a, which has been made plausi- 
ble on the grounds of a local entropy maximum (e.g. IH2007I) . 
However what has this to do with the black-hole bulge mass cor- 
relation that the co-eval growth was supposed to explain? 

If one assumes that the halo around the seed black hole is 
accreted eventually through two-body and clump-clump relax- 
ation ( MacMill an & Henr iksen 2003 ), then the mass of the black 
hole/halo should be proportional to the mass of the bulge simply 
because of the power law density (with nearly constant power in- 
side r s ). Letting M. be the black hole mass, M s be the bulge mass 
and r/,, r s , be the black hole halo and bulge radii respectively (r s 
is the scale radius); we find the relation 



M. 

M s 



(3- 2a) 



(29) 



If the appropriate a = 0.72, then the power of the radial depen- 
dence is 1.56. For a true universal relation however, we must 
assume that galaxy growth is not only self-similar, but that there 
is a similarity between galaxies so that the scale ratio is simi- 
lar. This is not exactly the case (e.g. Navarro et al., 20 09), but it 



2000 - 



1500 



1000 




-300 



-200 



-100 



Figure 3. The figure shows a contour plot of d 2 M/dEdJ for 
the full simulation of an isolated halo (MacMillan 2006). The 
separate plots of dM/dE and dMdJ are also indicated. The cor- 
responding density of states is almost independent of J and is the 

-2.5 



expected (e.g jBinney & Tremaine 1987" I \E\ law. The line 



on the upper left is J[E) for circular orbits. 

is nearly so. Such a relation might be tested by high resolution 
simulations containing a seed black hole. 

3.1. More general self-similar distribution functions 

The limits of equation ( fTBI ) that we discussed above are quite 
arbitrary, except for the expectation of central isotropy and the 
dominance of angular momentum in the outer regions. In this 
section we explore more general possibilities, even though em- 
pirically the preceding limits seem adequate. 

The distribution function ( fTBT ) can be put in a form that seems 
to generalize the FPDF (i.e. Fridmann and Polyachenko DF). We 

choose P oc k to obtain 



K 



U*)**\E„-E\W 



(30) 



The density integral over this DF requires a lower cut-off 
in angular momentum for a > 1 and an upper cut-off in energy 
for a < 1 in order to avoid singularities. We have exhibited the 
upper energy cut-off in (1301 explicitly. It appears as the arbitrary 
constant in the potential so that the difference E () — E still has 
the self-similar behaviour (r 2 ^ - ")), as does the density (r~ 2a ). 
Moreover for a > 1 the minimum angular momentum should be 
a fixed fraction of the value 2r 2 {E — <I>), say < k\ <C 1, in order 
for the density profile and potential associated with equation (TT3T > 
to apply. The radial velocity dispersion is also oc r2(i-«). it is 
only when a = 1 that one obtains the r~ 2 profile of the radial 
FPDF, and this must be discussed separately below. 

If one computes formally the density by integrating the DF 
( 1301 in a negative energy region (negative although a < 1, be- 
cause the central mass dominates), one finds that for a < 1 



(31) 
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As in the previous section, this can be used in an iterative fash- 
ion. 

Thus in a region dominated by a central point mass (and tak- 
ing E = 0) we see that the cusp would be, at lowest order, 



(32) 



The next order would be found by using this in the Poisson equa- 
tion to obtain a new potential and then a new density from equa- 
tion PTI l. Once again this low order result can not be flatter than 
r~ L5 .The steepest limit is r~ 2 as a — » 1—. We remark that at 
a = 2/3 one obtains the Bahcall and Wolf zero flux cusp r~ 7 / 4 , 
due to the filled loss cone in the DF. For a = 0.72, the cusp is like 
r~ 178 , very close to the Bahcall and Wolf cusp and still with a 
filled loss cone. In paper III of this series we shall find an exact 
example of this type. 

3.2. Singular isothermal sphere: the global inverse square law 

So far we have ignored the case when a = 1, the singular isother- 
mal limit. In fact the steady DF with a = 1 must be treated sep- 
arately because of the logarithmic potential <t> = V P In r (\P con- 
stant) that accompanies the inverse square density law. We may 
still deduce it from the characteristics of © if we take T = 
finally, since a steady state is the same at any time. The appro- 
priate characteristics become, with a = 1, 



— = 2P, 



dP 

ds 
dR 

ds 
dY 
ds 
dZ 
ds 



Y 

R, 

a 

1 / Z ¥ 



(33) 



— = 2Z. 



These integrate to give 

P = P(Z ,K)/Z, 
*P 

K = S InZ, 

2 



(34) 



where 



Li 

T 



z 

2P? 



^lnfi, 



and Z is again j 2 at s = on the characteristic. 

By writing the self-similar physical form at T = we find 
(for a = 1 ) 



P(K) 




j 2 










f 




2r 








J 2 




2r 



>Plnr In/ 2 

2 J 



(36) 



(37) 



Once again E = v 2 r /2 + j 2 /(2r 2 ) + y ¥\nr so that K = E - 
¥/21n j 2 . 

This limit with a = 1 is an extension of the inverse square 
density law to non-radial orbits, but the function of K is arbitrary. 
The DF is certainly non-unique for the same density profile. 



If we calculate the density profile for K > by integrating 
over ( f36b we arrive eventually at 



V2 



K 1 / 2 



dK 



ui 



du 



»i wyl — (u — V i'\nu)/2K 



(38) 



where u = j 2 /r 2 . The limits u\, 112 are the two roots found in 
terms of the Lambert function by setting the argument of the 
square root equal to zero. When fc/*P is large, the lower root 
becomes zero and the upper root approaches 2k. There is slow 
convergence of the inner integral. Provided that P(k) is such as 
to make the outer integral converge, the profile is r~ 2 . 

We note that for any finite K, the possible range of u has a 
definite upper limit and a definite lower limit that may be close 
to zero. This means that at given r there is a finite range in j 2 
that is present, and that this range collapses on zero as r — > 0. 
This is as expected with spherical symmetry. 

It is possible to have solutions with K < whereupon 



V2 



(>Pln 1 P- , P)/2 



P(k) d\K\ 
IkIV2 



"2 



du 



U\J ('PlnM 



■u)/2\k\-\ 
(39) 

where again the inner limits are the roots found by setting the 
argument of the square root equal to zero. For a real range of val- 
ues, fc| must be greater than a minimum value that depends 
on the value of l P. The same restriction on angular momentum 
as r — > applies. In all cases the velocity dispersion is constant. 

Unlike the radial FPDF that has an inverse square cusp, 
a black hole is not readily incorporated into this distribution. 
However at large enough r (or basically where in phase space 
where v r j r 2 <C v /r 2 ) it suffices to modify the integral by a point 
mass potential so that 



K = 



J 

2r 2 



■mar- 



GM. W 



2 J 



(40) 



The velocity dispersion is no longer strictly constant, unless 
j 2 ^> GM.r, but the influence of the black hole is weak in this 
regime. 

One special model of an inverse square, steady, anisotropic 
bulge is provided by equation ( l36l l by taking 



P(k) =Kexp 
(35) where b is any real number. 



"(f (1+*)) 



(41) 



In velocity space the anisotropic DF becomes explicitly 

,2 • 



(v z = v^ + v1) 



K 

— exp 



-( 



(b + l)\ 



) (vi)*- 



(42) 



This form is quite closely related to the DF suggested by Stiavelli 
and Bertin ( Stiavell i~& Bertin, 1985| l. However it offers a combi- 
nation of an exponential cut-off and a possible power-law cut-off 
in angular momentum. It does have the inverted distribution in 
energy (it increases outward) as suggested bv lMTJI This is also 
true of the DF that we suggest for the inner region, namely equa- 
tion ( fT8b with a < 1 . Thus we obtain this inversion naturally as 
a consequence of the self-similar evolution. 

When b = we have a pure Gaussian in the energy and once 
again the singular isothermal sphere. If b < there is a filled loss 
cone and for b > the loss cone is empty. Provided b > — 1 this 
DF is elliptical in velocity space (v r , v± ) at each r and becomes 
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circular (isotropic) at b = 0. The isothermal form with b = 
also appears in the fully asymmetric case of zeroth order coarse 
graining as discussed in the appendix of (H2004). 

This example is only of interest in the context of the persis- 
tent r~ 2 mass distribution in galaxies and clusters referred to in 
paper I There is clearly a lack of uniqueness however, as sev- 
eral distribution functions may produce an inverse square den- 
sity profile. Adding a black hole to this example is inconvenient 
(unless it is negligible) as the form of the DF is not suitable for 
iteration. However just as in the general case above, K may be 
taken as that given in equation d40t . 

One can calculate the growth of a central black hole embed- 
ded in an inner bulge with the DF ( TT8T >. We can expect to main- 
tain the steady bulge DF until the black hole mass is a significant 
fraction of the bulge mass. We find after obvious but tedious cal- 
culations that 

^=(2^/^^i,3 (1 - fl ) ) (43) 

at c z 

where x = r,/r s , the subscript s refers to characteristic bulge 
quantities, and a > 1/3 for convergence of the integrals. We re- 
call that the last stable radius of a Schwarzschild black hole is 
r. = 6GM,/c 2 . For a = 2/3 the growth is exponential, but the 
e-folding time is longer than the age of the Universe (5 x 10 19 s) 
even if the bulge has 1O 1O M in one kiloparsec. The result is 
similar for other values of a. However growth of halo cloaking 
the actual black hole could be much faster. In general the mass 
growth from the DF ([l"8l is (r/, and M/, are replaced by the black 
hole values to obtain the previous equation) 

^^Wp,^) 5 "-"', ,44, 

where M/, is the mass accreted inside the radius r/,. This radius 
would be an intermediate scale between the bulge scale r s and 
r., determined ultimately by the minimum angular momentum 
actually available. The accretion rate is proportional to the ra- 
dius of the structure and is faster than the black hole rate by 
r/j/r. when a = 1. One assumes that all of this mass is eventu- 
ally accreted by the black hole through relaxation processes (e.g. 
MacMilla n & Henri ksen 2003 ), although multiple steps might 
be required. Then the black hole growth is limited by the time 
to empty this reservoir onto the centre. 

This calculation is done without any bias towards the loss 
cone (/' 0), and with a self-similar DF of the form OOl l. one 
can actually grow the black hole directly. We shall save this cal- 
culation for paper HH1 

4. Conclusions 

In this paper, we have developed distribution functions that de- 
scribe both dark matter bulges and a central black hole or at least 
a central mass concentration. We succeed mainly in describing 
the dark matter bulges. 

In the discussion of cusps and bulges based on purely ra- 
dial orbits (paper|I]l, we were able to distinguish the Distribution 
function of Fridmann and Polyachenko from that of Henriksen 
and Widrow. The FPDF was found to describe accurately the 
purely radial simulations of isolated collisionless halos carried 
out in (MacMillan 2006). Moreover a point mass could be added 
without changing the form of the DF, which allows a self-similar 
growth of a central bulge or black hole. 

The final result concerning steady, self-similar radial 
orbits concerned the special case a = 1. The DF is a 



Gaussian that has been found previously in coarse graining 
dHenriksen & Le De lliou 2002). We included it there as a sec- 
ond example of a radial DF that produces an r~ 2 density profile 
(Mut ka 2009t . 

The inclusion of angular momentum here leads to more real- 
istic situations. We re-derived the steady self-similar DFs from 
first principles in equation d!5t . We showed that these can be 
used to describe the simulated collisionless halos calculated in 
(MacMill an 20061 1. if we use the value of a w 0.72 identified in 
(IH2007I) and take two limits. In one ( TT~8b the DF is isotropic and 
describes approximately the most relaxed central region of the 
bulge. The other limit dT9b describes the outer region and the 
transition across the NFW scale radius. The qualitative corre- 
spondence supports the relevancy of this family. The potential 
of a central mass can not be included exactly in these distribu- 
tion functions, but iteration is possible. An example was given 
for the DF COD. 

In particular, velocity dispersions have been calculated and 
compared qualitatively with the results of a high resolution nu- 
merical simulation of an isolated halo. These correspond to rea- 
sonable densities when adiabatic self-similarity is employed. 
Unfortunately these do not apply to the vicinity of the black 
hole itself. However the eventual dominance of the black hole 
through the r~ 1 potential becomes obvious, as seen in the iter- 
ated potential (I211l . 

Simple forms for the DF have been suggested for the various 
regions of the halo. And given the dynamic co-eval growth of 
black-hole/halo and bulge, we explain the black-hole bulge mass 
correlation simply d29l ). This relies on the approximate similarity 
between different halos. 

The forms that we find follow from the assumed self-similar 
path to equilibrium. Remarkably they are in accord with the 
general forms for the DF as presented in Stiavelli and Bertin 
flStiavelli & Bertin, 1985) and as modified in lMTJI That is, these 
distribution functions contain an inner inverted distribution in 
energy (referred to as negative temperature in MTJ) and a cut-off 
in the distribution in angular momentum farther from the centre. 
The scale of the inner region that is argued to coincide with NFW 
scale radius in this paper is the parameter r a in the earlier papers, 
so inner and outer have similar meanings. 

The variation is a power-law generally in both energy and 
angular momentum, although in the very relevant case of an in- 
verse square density profile, an exponential variation is possible. 
These conclusions follow from the assumed adiabatically self- 
similar evolution towards equilibrium. As such they tend to an- 
swer questions posed in lMTJI such as whether these distributions 
arise naturally and if the variations are exponential or power- 
law. We now know however that part of the answer lies in the 
anisotropy produced by the radial orbit instability dMWH 2 006). 

We have also found a self-similar generalization of the radial 
FPDF in equation ( l30t . Unfortunately this DF does not have the 
property of yielding a density that is independent of the potential 
and so a point mass potential is incompatible with self-similarity. 
It might describe a system at large radii with a large central mass. 

As always a = 1 must be treated separately and we give a 
derivation from first principles. The result is new and it gener- 
alizes the FPDF in the sense that p °< r~ 2 always. Although we 
can not place a central mass inside this bulge exactly, it allows 
an anisotropic DF in the r~ 2 bulge region. 

Generally we find that we can not describe elegantly 
anisotropic bulges containing black holes with self-similar DFs. 

In the next paper in this series (paperlHlb. we shall widen our 
scope to the study and production of non-self-similar cusps and 
DF. 
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